Régression spatiale

Author

Claude Grasland

INTRODUCTION

Objectifs

Cette séance va permettre :

  • de revoir plus en détail le concept de potentiel en utilisant le package R potential qui en facilite le calcul et la cartographie.

  • de revoir les modèles de régression multiple et de régression Poisson en les appliquant à la prédiction du nombre de clubs de sport présents ou absent d’une commune.

  • de mélanger des variables explicatives de type endogène (caractéristiques internes de la commune) et exogènes (situation de la commune par rapport aux espaces envrionnants)

Chargement des packages

L’installation est a priori la même que dans les sessions précédentes. On reprend juste pour mémoire la liste, à l’intention de ceux qui n’auraient pas suivi la séance précédente. On y ajoute un package pour le calcul de potentiel (potential) et un autre pour l’analyse des résultats de régression (car)

library(knitr)
library(dplyr)

library(sf)
library(mapsf)
library(RColorBrewer)
library(leaflet)
library(htmlwidgets)
library(htmltools)

library(ggplot2)
library(plotly)


library(potential)
library(car)

DONNEES

On reprend les trois jeux de données précédents en y ajoutant les deux bases de données sur le nombre de licences et le nombre de clubs présents dans chaque commune.

## Contour des communes
com<-readRDS("ParisPC/mapcom_parisPC.RDS")

## Population sur grille de 200 m
gri<-readRDS("ParisPC/gridpop_parisPC.RDS")

## Equipements localisés en latitude longitude
equ<-readRDS("ParisPC/equip_ParisPC.RDS") 

## Nombre de licences
lic <- readRDS("ParisPC/lic2018_ParisPC.RDS")

## Nombre de clubs
clu <-readRDS("ParisPC/clu2018_ParisPC.RDS")

Tableau de variables endogènes

On extrait de chaque fichier les informations relatives au sport retenu pour l’analyse et on centralise les résultats dans le fichier communal. On choisit de travailler sur l’exemple du football en prenant comme code de fédération 111 et comme code d’équipement 2802.

# Zone d'étude
tabcom <- st_drop_geometry(com)

# Extraction des clubs
myfede<-c("111")
myclu<- clu %>% filter(code_federation %in% myfede) %>%
                select(insee_com, nbclu =clubs_sportifs_2018)

# Extraction des licences du sport choisi
myfede<-c("111")
mylic<- lic %>% filter(code_fed %in% myfede) %>%
                group_by(insee_com) %>%
                summarise(nblic=sum(nblic))

# Extraction du nombre total de licences
lictot<- lic %>%  group_by(insee_com) %>%
                  summarise(lictot=sum(nblic))

# Extraction des équipements
typequ <- c("2802")
myequ <- equ %>% st_drop_geometry() %>%
                 filter(typ_code %in% typequ) %>%
                  group_by(insee_com) %>%
                  summarise(nbequ=n())

# Variables socio-eco
soceco<- gri %>% st_drop_geometry() %>%
                group_by(insee_com) %>%
                summarise(Ind = sum(Ind),
                          Men = sum(Men),
                          Men_pauv = sum(Men_pauv),
                          Ind_snv = sum(Ind_snv)) %>%
                mutate(pop = Ind,
                       men = Men,
                       txpauv = 100*Men_pauv/Men,
                       revhab = Ind_snv/Ind) %>%
                select(insee_com, pop,men, txpauv, revhab)
  

# Assemblage des tableaux
tabcom<- tabcom %>% left_join(soceco) %>%
                    mutate(denpop = pop/sup) %>%
                    left_join(lictot) %>%
                    left_join(mylic) %>%
                    mutate(pctlic = 100*nblic/lictot) %>%
                    left_join(myclu) %>%  
                    left_join(myequ) %>%
                     arrange(insee_com)

# Remplacement des NA par des 0
tabcom$nbclu[is.na(tabcom$nbclu)]<-0
tabcom$nblic[is.na(tabcom$nblic)]<-0
tabcom$nbequ[is.na(tabcom$nbequ)]<-0

# données + fonds de carte
mapcom <-left_join(com, tabcom)

CARTOGRAPHIE

Prenons l’exemple d’une cartographie du nombre de licenciés (stock) et de la part du total des licences (taux).

Avec mapsf

Normalement vous connaissez toutes les fonctions utilisez. Bien noter le changement de projection de la carte en crs = 2154.

mapcom<-mapcom %>% st_transform(2154)

mapdep<-mapcom %>% group_by(dept) %>% 
                summarise()

# Choix des classes 
    mybreaks<-c(0, 1,2,4,8,16,32,100)
# Choix de la palette 
    mypal<-brewer.pal(7, "RdYlBu")

# Carte de taux (choroplethe)
mf_map(mapcom, type="choro",
       var="pctlic",
       breaks= mybreaks,
       pal=mypal,
       border="white",
       lwd=0.5,
       leg_title = "% total licence",
       leg_val_rnd = 0,
       leg_pos = "topright")

# Ajout des départements
mf_map(mapdep, type="base",
       col=NA,
       border="black",
       lwd=1,
       add=T )

# Carte de stock (proportionelle)
mf_map(mapcom, typ = "prop",
       var = "nblic",
       inches = 0.05,
       col="gray50",
       leg_title = "nombre de licences",
       leg_pos = "topleft")

# Cadre, titre, ..
mf_layout(title = "Distribution des licences de golf dans Paris PC",
          credits = "Source : INSEE et Min. des Sports",
          frame = T,
          arrow=F,
          scale=T)

Avec leaflet

Programme plus compliqué … mais qui permet d’aboutir au même résultat avec une carte interactive. Noter que la projection n’est pas modifiée et demeure crs = 4326.

map<-mapcom %>% select(dept,insee_com, nom_com, nblic, pctlic ) %>% st_transform(4326)
map$lng<-st_coordinates(st_centroid(map))[,1]
map$lat<-st_coordinates(st_centroid(map))[,2]

mapdep<-mapcom %>% group_by(dept) %>% 
                summarise() %>%
                st_transform(4326)


# Choix de la variable
   myvar <-map$pctlic
# Choix des classes 
    mybreaks<-c(0, 1,2,4,8,16,32,100)
# Choix de la palette (c'est une fonction !)
   mypal <- colorBin('RdYlBu', 
                       myvar,
                       bins=mybreaks)
  
# Calcul du diamètre des cercles
   myradius <-8*sqrt(map$nblic/max(map$nblic, na.rm=T))  
   
# Préparation des popups
      mypopups <- lapply(seq(nrow(map)), function(i) {
      paste0(  paste("Commune               : ",map$nom_com[i]), '<br>',
               paste("Code INSEE           : " ,map$insee_com[i]), '<br>', 
               paste("Nb. de licences        : " ,map$nblic[i]), '<br>', 
               paste("% total licence    :", round(map$pctlic[i],1))
            ) 
            })
      mypopups<-lapply(mypopups, htmltools::HTML)





map <- leaflet() %>% 
            addProviderTiles('Esri.WorldTopoMap') %>%

  # Réalisation de la carte choroplèthe
            addPolygons(data = map,
                        fillColor = ~mypal(pctlic),
                        fillOpacity = 0.5,
                        color = "white",
                        popup = mypopups,
                        weight = 1,
                        highlightOptions = highlightOptions(weight = 3, color = 'green')) %>%

  # Ajout de la carte des départements
              addPolygons(data = mapdep,
                        fill = FALSE,
                        color = "black",
                        weight = 2) %>%

  # Ajout de la carte de stocks    
               addCircleMarkers(data=map,
                              lat = ~lat,
                              lng = ~lng,
                              radius = myradius,
                              stroke = FALSE,
                              label = ~nblic,
                              fillColor = "gray50",
                              fillOpacity = 0.5)%>%
   
  # Ajout de la légende 
            addLegend(data = map,
                      pal = mypal, 
                      title = "% licences",
                      values = ~pctlic, 
                      position = 'topright') 



map

REGRESSION LINEAIRE

On va ici un modèle de régression n’utilisant que des variables endogènes (internes aux communes). On essaye de modéliser Y = % de licences pour le sport considéré

Choix des variables

Au vu de la carte précédente, on peut penser à plusieurs variables explicatives et faire des hypothèses sur le résultat attendu

  • X1 : densité de population (relation positive car il faut de la place pour installer un terrain de football)
  • X2 : revenu moyen par habitant (relation négative car les riches pratiquent plutôt des sports d’élites ce qui n’est pas le cas du football)
  • X3 : % de ménages pauvres (relation positive car le football est considéré à tort ou à raison comme un outil de promotion sociale)

On crée un tableau avec Y, X1, X2, X3

don<-tabcom %>% select(insee_com, nom_com, Y=pctlic, X1=denpop, X2=revhab,X3=txpauv)
head(don)
  insee_com                  nom_com        Y        X1       X2       X3
1     75101 PARIS-1ER-ARRONDISSEMENT 1.703085  8391.209 38433.52 13.23255
2     75102  PARIS-2E-ARRONDISSEMENT 5.564024 21350.505 34985.14 15.00129
3     75103  PARIS-3E-ARRONDISSEMENT 3.400690 30778.632 36195.03 13.99706
4     75104  PARIS-4E-ARRONDISSEMENT 4.714971 15908.750 36304.64 13.55867
5     75105  PARIS-5E-ARRONDISSEMENT 3.075031 19597.619 37720.86 13.18044
6     75106  PARIS-6E-ARRONDISSEMENT 4.521818 16633.721 43607.03 12.49411

Analyse des corrélations

On analyse la matrice de corrélation :

cor(don[,3:6])
            Y          X1         X2          X3
Y   1.0000000 -0.35258137 -0.8156021  0.72652758
X1 -0.3525814  1.00000000  0.2526739  0.04524385
X2 -0.8156021  0.25267387  1.0000000 -0.72002742
X3  0.7265276  0.04524385 -0.7200274  1.00000000

On vérifie la forme des relations et on teste leur significativité

  • Densité de population
scatterplot(don$X1,don$Y)

cor.test(don$X1,don$Y)

    Pearson's product-moment correlation

data:  don$X1 and don$Y
t = -4.474, df = 141, p-value = 0.00001569
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.4884600 -0.2000084
sample estimates:
       cor 
-0.3525814 
  • Richesse par habitant
scatterplot(don$X2,don$Y)

cor.test(don$X2,don$Y)

    Pearson's product-moment correlation

data:  don$X2 and don$Y
t = -16.738, df = 141, p-value < 0.00000000000000022
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.8640692 -0.7521516
sample estimates:
       cor 
-0.8156021 
  • % ménages pauvres
scatterplot(don$X3,don$Y)

cor.test(don$X3,don$Y)

    Pearson's product-moment correlation

data:  don$X3 and don$Y
t = 12.555, df = 141, p-value < 0.00000000000000022
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.6385290 0.7957734
sample estimates:
      cor 
0.7265276 

Toutes les variables affichent des corrélations linéaires significatives ! Mais on note que la forme des relations n’est pas forcément linéaire. On y reviendra …

Modélé additif

On teste le modèle additif suivant

\(Y = a_0+a_1.X_1+a_2.X_2+a_3.X_3+ \epsilon\)

modlin<-lm(Y ~ X1 + X2 + X3, data=don)
summary(modlin)

Call:
lm(formula = Y ~ X1 + X2 + X3, data = don)

Residuals:
    Min      1Q  Median      3Q     Max 
-7.7420 -1.9405 -0.2735  1.5522 12.0714 

Coefficients:
               Estimate  Std. Error t value             Pr(>|t|)    
(Intercept) 19.59836228  2.06344477   9.498 < 0.0000000000000002 ***
X1          -0.00020419  0.00003642  -5.607       0.000000107254 ***
X2          -0.00039364  0.00005677  -6.934       0.000000000141 ***
X3           0.35444421  0.05503644   6.440       0.000000001810 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.117 on 139 degrees of freedom
Multiple R-squared:  0.7598,    Adjusted R-squared:  0.7546 
F-statistic: 146.6 on 3 and 139 DF,  p-value: < 0.00000000000000022
anova(modlin)
Analysis of Variance Table

Response: Y
           Df Sum Sq Mean Sq F value                Pr(>F)    
X1          1  699.1   699.1  71.940   0.00000000000002938 ***
X2          1 3170.6  3170.6 326.279 < 0.00000000000000022 ***
X3          1  403.0   403.0  41.476   0.00000000180962499 ***
Residuals 139 1350.7     9.7                                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  • Conclusion : toutes les variables sont significatives et affichent les signes attendus. L’analyse de variance montre cependant que c’est la variable X2 (revenu par habitant) qui est la plus déterminante. La qualité de l’ajustement est élevée (r2 = 76%)

Modèle multiplicatif

Essayons maintenant un modèle multiplicatif de la forme

\(log(Y) = a_0+a_1.X_1+a_2.X_2+a_3.X_3+ \epsilon\)

qui correspond à une forme multiplicative de l’effet des variables puisque l’on a :

\(Y = exp(a_0+a_1.X_1+a_2.X_2+a_3.X_3+ \epsilon)\)

et donc :

\(Y = exp(a_0).exp(a_1.X_1).exp(a_2.X_2).exp(a_3.X_3)\)

modexp<-lm(log(Y) ~ X1 + X2 + X3, data=don)
summary(modexp)

Call:
lm(formula = log(Y) ~ X1 + X2 + X3, data = don)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.2206 -0.1266  0.0045  0.1464  0.5606 

Coefficients:
                Estimate   Std. Error t value             Pr(>|t|)    
(Intercept)  3.817787144  0.157533985  24.235 < 0.0000000000000002 ***
X1          -0.000012995  0.000002781  -4.673           0.00000692 ***
X2          -0.000053438  0.000004334 -12.330 < 0.0000000000000002 ***
X3           0.007414592  0.004201765   1.765               0.0798 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.238 on 139 degrees of freedom
Multiple R-squared:  0.7949,    Adjusted R-squared:  0.7905 
F-statistic: 179.6 on 3 and 139 DF,  p-value: < 0.00000000000000022
  • Conclusion : Pas de différences dans le signe des coefficients et leur significativité. Mais la qualité de l’ajustement est plus élevée (r2 = 79%) ce qui signifie que le modèle multiplicatif décrit plus fidèlement les effets constatés.

POTENTIELS

On va maintenant utiliser le package potential pour créer des variables de type exogène mesurant la distribution des indicateurs non pas dans la commune mais dans le voisinage de celle-ci.

Potentiel d’équipement

# Extraction des équipements et projection 2154
typequ <- c("2802")
myequ <- equ %>%filter(typ_code %in% typequ) %>%
                  mutate(nb = 1) %>%
                  st_transform(2154)

# projection des communes en 2154
mapcom <- mapcom %>% st_transform(2154)

# Distance communes-équipement
dist<-create_matrix(myequ,mapcom)

# calcul du potentiel
mapcom$pot_equ_2000<-potential(x=myequ,   # Ressources
                               y=mapcom,  # Population
                               d=dist,    # Distance population x Ressources
                               var = "nb", # Quantité de ressource
                               fun="e",    # famille exponentielle
                               span = 2000, # distance ou f(Dij) = 0.5
                               beta=2)      # Exposant de la distance


# Cartographie du résultat
mf_map(mapcom, type="choro",var="pot_equ_2000")
mf_map(myequ, type="base",add=T,col="red")

  • Commentaire : La carte montre que les arrondissements du centre de Paris qui ne disposent pas de terrain de football ont cependant un potentiel d’acès à ceux-ci dans un voisinage gaussien de 2 km.

Potentiel de licences

# Extraction des licences projection 2154
maplic<-mapcom%>% select(nblic) %>%st_transform(2154)

# projection des communes en 2154
mapcom <- mapcom %>% st_transform(2154)

# Distance communes-équipement
dist<-create_matrix(maplic,mapcom)

# calcul du potentiel
mapcom$pot_lic_2000<-potential(x=maplic,y=mapcom,d=dist,var = "nblic", fun="e",span = 2000,beta=2)

# Cartographie du résultat
mf_map(mapcom, type="choro",var="pot_lic_2000")
mf_map(maplic, type="prop",var="nblic",col="red", inches=0.05)

  • Commentaire : La carte montre que le potentiel de sportifs ayant une licence de football est maximale dans les arrondissements périphériques de Paris et dans la banlieue Nord.

REGRESSION DE POISSON

On va construire un modèle de régression de Poisson pour prévoir le potentiel de clubs du potentiel d’équipement et de son potentiel de licence.

Choix des variables

Le tableau comporte uniquement des stocks et des potentiels

Y : nombre de clubs X1a : nombre d’équipements dans la commune X1b : potentiel d’équipement dans un voisinage gaussien de 2.5km X2a : nombre de joueurs ayant une licence dans la commune X2b : potentiel de joueurs ayant une licence dans un voisinage gaussien de 2.5km

On crée un tableau avec Y, X1a, X1b X2a,X2b

don<-mapcom %>% select(insee_com, nom_com, 
                       Y=nbclu, 
                       X1a = nbequ,
                       X1b = pot_equ_2000, 
                       X2a = nblic,
                       X2b = pot_lic_2000) %>%
  st_drop_geometry()
head(don)
  insee_com                  nom_com  Y X1a       X1b  X2a       X2b
1     75119 PARIS-19E-ARRONDISSEMENT 17   3 37.262588 2162 22122.189
2     75106  PARIS-6E-ARRONDISSEMENT  5   0  6.116627  400 14617.246
3     93071                   SEVRAN  3  13 47.018154 1190 10356.941
4     94056                  PERIGNY  2   1  5.746092   84   814.633
5     94041           IVRY-SUR-SEINE  5  10 38.133234 1305 17336.098
6     75112 PARIS-12E-ARRONDISSEMENT 14  29 65.061731 1393 28111.675

Analyse des corrélations

On analyse la matrice de corrélation :

cor(don[,3:7])
            Y       X1a       X1b       X2a       X2b
Y   1.0000000 0.4315624 0.2416476 0.7516748 0.5794876
X1a 0.4315624 1.0000000 0.6279578 0.5791498 0.2387389
X1b 0.2416476 0.6279578 1.0000000 0.5980269 0.4025822
X2a 0.7516748 0.5791498 0.5980269 1.0000000 0.5992941
X2b 0.5794876 0.2387389 0.4025822 0.5992941 1.0000000

On vérifie la forme des relations et on teste leur significativité

  • Equipement de la commune
scatterplot(don$X1a,don$Y)

cor.test(don$X1a,don$Y)

    Pearson's product-moment correlation

data:  don$X1a and don$Y
t = 5.6808, df = 141, p-value = 0.00000007394
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.2878018 0.5563023
sample estimates:
      cor 
0.4315624 
  • Potentiel d’équipement dans le voisinage de la commune
scatterplot(don$X1b,don$Y)

cor.test(don$X1b,don$Y)

    Pearson's product-moment correlation

data:  don$X1b and don$Y
t = 2.957, df = 141, p-value = 0.003643
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.08070004 0.39031390
sample estimates:
      cor 
0.2416476 
  • Nombre de licences dans la commune
scatterplot(don$X2a,don$Y)

cor.test(don$X2a,don$Y)

    Pearson's product-moment correlation

data:  don$X2a and don$Y
t = 13.533, df = 141, p-value < 0.00000000000000022
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.6702226 0.8152346
sample estimates:
      cor 
0.7516748 
  • Potentiel de licences dans le voisinage de la commune
scatterplot(don$X2b,don$Y)

cor.test(don$X2b,don$Y)

    Pearson's product-moment correlation

data:  don$X2b and don$Y
t = 8.4432, df = 141, p-value = 0.00000000000003385
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.4590001 0.6790441
sample estimates:
      cor 
0.5794876 

Modèle

On teste le modèle suivant :

\(Y = exp(\alpha+\beta_{1a}.X_{1a}+\beta_{1b}.X_{1b}+\beta_{2a}.X_{2a}+\beta_{2b}.X_{2b} +\epsilon)\)

On emploie une régression de Poisson car le nombre de clubs est une variable quantitative discrète composée d’entier. On suppose qu’elle obéit à une loi de Poisson et on vérifie si c’est exact en comparant la moyenne et l’écart-type de Y.

mean(don$Y)
[1] 3.79021
sd(don$Y)
[1] 3.711143

Des tests plus précis pourraient être utilisés mais dans le cas présent on voit que les deux valeurs sont très proche et qu’il est acceptable de considérer que le nombre de club obéit bien à une loi de Poisson.

modpoi<-glm(Y ~ X1a + X1b + X2a + X2b, data=don, family="poisson")
summary(modpoi)

Call:
glm(formula = Y ~ X1a + X1b + X2a + X2b, family = "poisson", 
    data = don)

Coefficients:
               Estimate  Std. Error z value             Pr(>|z|)    
(Intercept)  0.25253750  0.16302149   1.549             0.121356    
X1a          0.03330562  0.00952547   3.496             0.000471 ***
X1b         -0.02185387  0.00406491  -5.376         0.0000000761 ***
X2a          0.00086511  0.00008965   9.650 < 0.0000000000000002 ***
X2b          0.00005579  0.00001095   5.096         0.0000003464 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 417.02  on 142  degrees of freedom
Residual deviance: 121.44  on 138  degrees of freedom
AIC: 534.46

Number of Fisher Scoring iterations: 5
Anova(modpoi,type="III")
Analysis of Deviance Table (Type III tests)

Response: Y
    LR Chisq Df            Pr(>Chisq)    
X1a   11.676  1             0.0006331 ***
X1b   29.632  1         0.00000005223 ***
X2a   94.633  1 < 0.00000000000000022 ***
X2b   25.540  1         0.00000043329 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Discussion

Les quatre facteurs conditionnent bien l’apparition des clubs de football dans les communes et sont tous significatifs. Le nombre de terrain de football dans la commune a un effet positif sur l’apparition d’un ou plusieurs club dans la commune ce qui parait logique. Par contre la quantité de terrains de football dans les communes voisines a un effet négatif ce qui peut se comprendre comme un effet de concurrence. Quant au nombre de personnes détenant une licence, il a un effet positif à la fois dans la commune et dans les communes voisines.

La régression de Poisson n’offre pas directement de qualité d’ajustement comme dans le cas de la régression par la méthode des moindres carrés ordinaires. On peut néanmoins calculer un pseudo-R2 dit de Mc Fadden en calculant la part de déviance expliquée c’est-à-dire en effectuant le calcul suivant :

\(R^2_{McFadden} = \frac{NullDeviance-ResidualDeviance}{NullDeviance}\)

(417.02-121.44)/(417.02)
[1] 0.7087909

Dans notre exemple on trouve un \(R^2\) de Mc Fadden d’environ 71% ce qui est très satisfaisant et montre que la connaissance de la localisation des terrains et des licences permet dans une large mesure de prédire le nombre de clubs présents dans une commune.

PROLONGEMENTS

Pour en savoir plus sur les modèles de régression, le plus simple est de partir du petit billet introduction au GLM de Claire Della Vedova qui explique bien la différence entre le modèle linéaire classique et les modèles de régression logistique ou de régression de Poisson.

Ensuite … il faudra s’attaquer à un bon manuel de statistiques. Nous recommandons celui de Daniel J. Denis Univariate, Bivariate Multivariate Statistics using R qui est très pédagogique et fournit les programmes d’application en R ou en Python.